/*!  \file  BiReduction.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      09/19/2013 02:19:16 PM
 *   \copyright GNU Public License.
 */
#include <float.h>
#include <stdio.h>
#include "UnitTester.h"
#include "Matrix.h"
#include "ImplementationFactory.h"
#include "BiReduction.h"

BiReduction::BiReduction(const char * name_): Implementation(name_, BiReductionGroup){
	setCompileTime(__DATE__, __TIME__);
};
BiReductionUnitTest::BiReductionUnitTest(BiReduction & reduct_): reduct(reduct_) {
	int l = strlen(reduct.getName());
	char buf[l + 20]; 
	sprintf(buf, "BiReduction %s", reduct.getName());
	setName 				(buf);
	setDescription  	("Test registered bi reduct implementations");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	setCompileTime(__DATE__, __TIME__); // each class that inherits this one should also have this line.
	setCompileTime(reduct.getCompileTime());
};
bool  BiReductionUnitTest::test1(const int m, const int n, bool transA, const double limit) const {
	assert(m * n <= 30);	
	double a[] = {8.79,  6.11, -9.15,  9.57, -3.49,  9.84,
	              9.93,  6.91, -7.93,  1.64,  4.02,  0.15,
	              9.83,  5.04,  4.86,  8.83,  9.80, -8.99,
	              5.45, -0.27,  4.85,  0.74, 10.00, -6.02,
	              3.16,  7.98,  3.01,  5.80,  4.27, -5.31};
	Matrix A   (n, m, General, a); 
	if (transA) 
		A.transpose();
	Matrix tmp (A, CopyWithData);

	Matrix B   (A.rows(), A.cols(), BiDiagonal );
	Matrix Bgen(A.rows(), A.cols(), General); Bgen.fill(0); 
	Matrix Q   (A.rows(), A.rows(), ElementaryReflectors);
	Matrix P   (A.cols(), A.cols(), ElementaryReflectors);

	reduct.compute(tmp, Q, B, P ); // tmp is destroyed in this computation
	Bgen += B;
	const int nNorms = 4;
	double norm[4];
	double normA = A.nrm2();
	norm[0] = (        A     - Q * Bgen * P.T()).nrm2();
	norm[1] = (Q.T() * A     -     Bgen * P.T()).nrm2();
	norm[2] = (        A * P - Q * Bgen        ).nrm2();
	norm[3] = (Q.T() * A * P -     Bgen        ).nrm2();
	bool res = true;
	for (int i = 0; i < 4; i++) {
		res &= norm[i] < max(n,m) * normA * DBL_EPSILON * limit;
	}
	return res;
	

}
bool BiReductionUnitTest::test () const {
	bool res = true;
	double limit = 1.0;
	res &= test1(5, 5,  true, limit);
	res &= test1(6, 5,  true, limit);
	res &= test1(5, 6,  true, limit);
	res &= test1(5, 5, false, limit);
	res &= test1(5, 6, false, limit);
	res &= test1(6, 5, false, limit);
	res &= test1(3, 6,  true, limit);
	res &= test1(6, 4, false, limit);

	return res;
}
